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Abstract 

Due to discretization effects and truncation to finite domains, many electromagnetic simulations present 
non-physical modifications of Maxwell’s equations in space that may generate spurious signals affecting 
the overall accuracy of the result. Such modifications for instance occur when Perfectly Matched Layers 
(PMLs) are used at simulation domain boundaries to simulate open media. Another example is the use 
of arbitrary order Maxwell solver with domain decomposition technique that may under some condition 
involve stencil truncations at subdomain boundaries, resulting in small spurious errors that do eventually 
build up. In each case, a careful evaluation of the characteristics and magnitude of the errors resulting 
from these approximations, and their impact at any frequency and angle, requires detailed analytical and 
numerical studies. To this end, we present a general analytical approach that enables the evaluation of 
numerical discretization errors of fully three-dimensional arbitrary order finite-difference Maxwell solver, 
with arbitrary modification of the local stencil in the simulation domain. The analytical model is validated 
against simulations of domain decomposition technique and PMLs, when these are used with very high-order 
Maxwell solver, as well as in the infinite order limit of pseudo-spectral solvers. Results confirm that the new 
analytical approach enables exact predictions in each case. It also confirms that the domain decomposition 
technique can be used with very high-order Maxwell solver and a reasonably low number of guard cells with 
negligible effects on the whole accuracy of the simulation. 
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1. Introduction 


Very high-order Finite-Difference Time-Domain (FDTD) Maxwell solver and, in the infinite order limit 
pseudo-spectral solvers, are methods of choice for solving three-dimensional electromagnetic problems that 
require very high accuracy over a large band of frequencies. In Particle-In-Cell (PIC) simulations of laser- 
plasma mirror interactions for instance, a good dispersion relation is needed over a large band of frequencies 
and angles for the modeling of relativistic harmonic generation mm- In PIC simulations of laser-plasma 
acceleration, numerical dispersion induces numerical Cherenkov effects mm that can be highly detrimental 
for the modeling of ultra relativistic beams and plasmas. In general, the mitigation of all these effects 
requires the use of very high spatial and temporal resolutions which, with current low order FDTD solvers, 
sometimes prevent parametric studies of large parameter space at scale in full three dimensions. In the case 
of numerical Cherenkov effects, it was shown analytically and numerically that pseudo-spectral solvers are 
generally more stable than standard second-order schemes j5]. The use of very high-order or pseudo-spectral 
solvers can significantly decrease the resolution needs and increase the overall stability for a given accuracy 
and can thus enable realistic 3D PIC simulation studies that are otherwise not practical. 

Despite significant advantages in terms of accuracy and memory-minimization |^, high-order and pseudo- 
spectral solvers have however not been widely adopted so far for large-scale simulations on massively parallel 
supercomputers because of their low parallel scalability, which is a direct consequence of their spatial non¬ 
locality (large stencils). Indeed, these solvers commonly use Fast Fourier Transform (FFT)-based algorithms 
that require global inter-processor communications in the computation of Fourier transforms, limiting their 
scaling to a few thousands of cores. This is preventing the use of pseudo-spectral solvers in very computation¬ 
ally demanding 3D plasma electromagnetic simulations that usually mandate several hundred of thousands 
of cores. 

Recently, a new method jTj for solving time-dependent problems (e.g. Maxwell’s wave equations) proposed 
to apply the cartesian domain decomposition technique currently used with low order FDTD solvers to 
very high-order pseudo-spectral solvers. As in the case of low-order schemes, the simulation domain is 
divided into several subdomains and Maxwell’s equations solved locally on each subdomain using local 
convolution (e.g. FDTD), or local FFT’s when the order p becomes very large. The fundamental argument 
legitimating this method is that physical information cannot travel faster than the speed of light. Choosing 
large enough guard regions should therefore ensure that spurious signal stemming from stencil truncations at 
subdomains boundaries is practically limited to the guard regions and that only a negligible fraction reenters 
the simulation domain after one time step. One potential drawback of this approach is the non-locality of 
Gibbs oscillations arising when a signal is truncated at the edges of the guard regions. Indeed, depending 
on the size of guard regions, errors stemming from stencil truncations may affect the entire simulation 
domain and instabilities could build up. Studying the impact of stencil modifications/truncations on the 
overall accuracy of the simulation is thus crucial to validate the use of arbitrary order schemes with domain 
decomposition techniques. 

1.1. Goals and outline of this study 

The goal of this study is to provide a very general analytical approach that enables the prediction of 
the total amount of error coming from stencil modifications in 3D electromagnetic simulations. For a fixed 
simulation configuration, this approach can be used to compute the optimal choice of numerical parameters 
(e.g. solver order p, space and time steps, number of guard cells) that will compute the solution in a 
minimum time and with a guaranteed accuracy. Our model will be especially relevant to predict errors when 
cartesian domain decomposition is used along with very high order/pseudo-spectral solvers in electromagnetic 
simulations. It is also applicable to high-accuracy prediction of the performance of PMLs with high-order or 
pseudo-spectral solvers. 

The study is divided in three sections: 

(i) Section 1 presents a very general method that analytically calculates the error induced by any modifi¬ 
cations of Maxwell’s equations in the simulation domain. 
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(ii) In sections 2 and 3, our analytical model is benchmarked against simulations. In particular, it is shown 
that our model can be used to deduce the growth rate of spurious truncation signals induced by the 
domain decomposition (section 2) as a function of various numerical parameters (Number of guard 
cells, order of Maxwell solver, subdomain sizes). It is also demonstrated that the new model enables 
very accurate predictions that are superior to previous models, of the efficiency of PMLs (section 3) 
with high-order stencils as a function of frequency and angle. 

2. Model for errors induced by spatial modifications of Maxwell’s equations in the simulation 
domain 

This section presents a method for deriving analytically the error induced by the modifications of 
Maxwell’s equations at an arbitrary number of nodes of the grid in the case of a plane and monochro¬ 
matic wave. The total error for an arbitrary waveform/beam is derived using plane wave decomposition and 
linearity of Maxwell’s equations. 

When the simulation domain is uniform, the Von Neumann analysis [8] provides a very simple yet powerful 
means of studying the stability and growth of errors in the simulation domain due to the discretization of 
equations. However, as soon as the simulation domain has discontinuities at any point (e.g at domain 
boundaries, mesh irregularities or stencil variations) where the discretized Maxwell’s equations are modified, 
it is no longer adapted. 

In the following we provide an alternative approach that provides accurate estimates of the total error 
( induced by the modifications of Maxwell’s equations at arbitrary nodes of the simulation domain. In this 
paper, we consider the most common case of a Maxwell’s field solver on a staggered grid [^. In appendix 
A, we give a formal definition of this scheme at order p and verify that when p —)• oo, this solver converges 
to the pseudo-spectral solver |10| . verifying that the new analysis model applies to pseudo-spectral solvers 
when p —)■ oo. 

5.1. Principle of the technique 

The principle of our analytical approach is progressively introduced in three steps by considering solutions 
of Maxwell’s equations for the three different cases: 

(i) regular stencil in vacuum, 

(ii) regular stencil with an external monochromatic source point, 

(iii) irregular stencil in vacuum (case of interest). 

2.1.1. Regular stencil in vacuum 

Discretized Maxwell’s equations in vacuum are written as: 

= E" - c6tV* X (1) 

-cJtV; X cE" (2) 

where n is the time step, E the electric field and B the magnetic field defined on staggered grids r and 
s. V* is the discrete finite difference operator of order p for the staggered scheme. 

A plane wave of frequency lv will propagate with a wavevector k given by the dispersion relation k = /(cu), 
where k is the norm of k. In ID, this relation is: 

sin = px Cf sm{2l — l)k— 

1=1 

where 5t and 5x are respectively the time step and mesh size, the coefficients Cf are the stencil coefficients 
of the order p staggered scheme solver (cf. Appendix A) and px = c5t/5x is the Courant parameter. 


( 3 ) 
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2.1.2. Regular stencil with an external monochromatic source point term 

We now investigate the discrete solutions of Maxwell’s equations in ID when an external monochromatic 
source point J”" = with = —1, is introduced at position i = 0 on the x-axis: 


n+1 _ 




= E” - c6tV*p X cB, 


n+t 

«+2 ^ jn 


= Bs 


" - c6tV*p X cE" 


( 4 ) 

( 5 ) 


As Maxwell’s equations are linear, all the generated electromagnetic fields will also be monochromatic 
with frequency a; and Maxwell’s equations (4j5) can be written: 


g. ( gici^W2 _ 


& _i_i 


p/2 

1=1 

p/2 

-rjx E CP (e,+, - 
9=1 


( 6 ) 

( 7 ) 


where we used EP = and . Replacing the expressions of and 

given by equation Q in equation yields the following equation for the electric field: 


ui6t uj5t 

4sin2 + ??^ ^ Cf Cf [ei+(g+i)_i - (ei+(,j_;) + ei_(^q_i)) + ei_(q+,)+i] = -2jc sin —Co5(/) (8) 

d,9)=i 

which can be written in a more compact form: 


p-i 


Lj6t 


^ Kf ei+i = -2jc sin — 


(9) 


i=-p+i 


where the coefficients are functions of the stencil coefficients rjx and uiSt. These coefficients are 
symmetric and verify Kf = kF_,. 


Calculating the electric field Cj at node i thus requires the values of the electric field at 2p — 1 adjacent 
nodes, as illustrated by the black dotted line on Fig. The system of equations written at all nodes i of 
the domain is closed by considering that at positions /+ > 0 and i_ < 0 far from the external source point, 
equation ([^ corresponds to a propagation equation in vacuum (cf. Fig[^ (a) and (c)) and admits solutions 
of the form: 


^i± — 


_ ^^^±jci±k5x 


( 10 ) 


where k is solution of the dispersion relation k = f{uj) of the scheme in vacuum]^ and the amplitude 
of the radiated wave in vacuum. Unrolling the system of equations is equivalent to "retropropagating" 
the single waves given by equation (10) from the far field to the initial source point i = 0. 

For order p = 2, equation ([^ can be solved independently on subdomains (i < 0) and (i > 0) (cf. Fig[^ 
(a) and (b)) and corresponds to a simple propagation in vacuum from both sides of the source point: 


r i ^ 0 

\ i ^ 0 


( 11 ) 
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Figure 1: Field emitted by a monochromatic external source point of amplitude ^o- In each panel, the black dotted line represents 
the envelope of the stencil with coefficients Ki used to compute d from equation § at node i. The external monochromatic 
source is represented in red. Adjacent secondary sources, created by stencils crossing on each side of the primary source location 
for p > 2, are represented in green and blue, (a) Calculation of the field a far from the external source point for order p — 2. 
(b) Calculation of the field near the external source point for order p = 2 (c) Calculation of the field ei far from the external 
source point for order p > 2. (d) Calculation of held d near the external source point for order p > 2. 


with: 


— 




2rix cos ^ 


( 12 ) 


For order p > 2, the solution far from the source point (cf. Fig. [^(c)) is still given by equation (10) and 
corresponds to a plane monochromatic wave propagating with wavevector k given by the dispersion relation 
of the order p scheme. 

However, near the external source point (i.e —{p — 1) < i < p — 1), equation ([^ does not correspond 
to a propagation equation anymore and its solutions will not be the same as for order 2, i.e single waves 
propagating in the backward (for i<0) or forward (for i>0) directions. 

Indeed, equations on the half domains i < 0 and i > 0 are now coupled due to a larger spatial extent 
of the high order stencil (cf Fig. (d)) and this coupling can be interpreted as the presence of secondary 
source points near the initial external source point. 

Considering the symmetry of the external current distribution with respect to i = 0, the solutions e* = e_i 
are also symmetric and equation ([^ written at positions —(p — 1) < io < (p — 1) (ig 0) near the source 
point thus can be written: 


(p-i) 

^ Kf*ei^+i = 0 

/=-(p-i) 


(13) 


where Kf* is different from its value Kf in vacuum. For the backward case — (p — 1) < io < 0, Kf* is: 


— (p — 1) < i + I < —i — (p — 1) 


p* _ 


= <{ 2Kf -i-(p-l) <i + l<0 (14) 

0<i+l<i+p—1 

While the solution given by equation © satisfies the propagation equation ([^ in vacuum for ig ^ 
— (p — 1) or io ^ (p — 1), it does not satisfy anymore equation (13) as soon as io ^ —(p — 2) or ig ^ (p — 2). 


This modification of the electric field equation will thus change the value of the field at io compared to 


a simple propagation in vacuum as given by equation (11). This can be formally written as adding another 


monochromatic "secondary source" point in i = io of amplitude that will modify the amplitude of the 
solution in vacuum. 


5 



























As for the original external source point, this secondary source point also radiates an electromagnetic 
wave in i < 0 and i > 0 directions and thus contribute to the total electromagnetic held far from its position 
(cf Fig. (d)). The total electric held thus has the following form: 


e,- oc 




^i+l<p 


i ^i- = -{p - 1) 
i^i+ = (p-l) 


(15) 


Far from the external source, the total electric held given by equation (15) can still be written as a 


single wave propagating forward (for i ^ i_|_) or backward (for i ^ i_) . This single wave results from the 
interference of multiple waves generated by the initial source point and by the 2{p — 2) surrounding secondary 


sources. 


However, near the external source { — {p — 2) ^ i ^ {p — 2)), the form of the total electric held ( |15[ ) is 
more complex and is a linear combination of multiple waves propagating forwards and backwards (cf Fig. 

(d)). 


2.1.3. Irregular stencil in vacuum 

In the following, we show that modifying Maxwell’s equations at one node is equivalent to adding a 
"pseudo-source" term. In accordance with the result of the external source case described in the previous 
subsection, this induces the creation of additional adjacent secondary pseudo-sources around the modihed 
node. 

Formally, local modihcations of the stencil can be written as: 


E^+i 


- l3rc5tV** X B"^^ 

(16) 

n-i 

a,Bs 2 - l3scdtV*p* X cE” 

(17) 


where a, /3 and V** are space varying coefficients and operator. 

Modifying the electric held equation at one node ro = i'f'o,x,'i"o^y,ro^z) only of the simulation domain can 
thus be written as: 

En+i = E" - c6tV*p X + J” (18) 

which is the same equation as in vacuum but with an additional time-dependent pseudo-source term J**: 


J” = 


(-1 + a,)E(l - /3rc5tiV** - V!) X B 


n+^ 


6(r - ro) 


(19) 


For a monochromatic driving wave of frequency ui propagating with wavevector k given by the dispersion 
relation in vacuum (far from the modihcation), the modihcation of Maxwell’s equations at rg will thus 
generate a monochromatic pseudo-current J** at r = tq (cf. Fig. I (a)): 

J" = - ro) (20) 


Similarly to the previous case of an external point source term at order p > 2, this "pseudo-source" term will 
induce the creation of {2p — 2) monochromatic secondary source terms of amplitudes (Cro+l)-{p- 2 )^l^{p- 2 ) 
each dimension around the modihed node rg, that will also radiate and contribute to the total electric held. 
These coefficients will be further called re-emission coefficients. 

A similar reasoning on B shows that 2{p — 2) time-dependent secondary pseudo-source terms are produced 
around modihed internodes. 
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2.1.4- New features added by the model 

Former techniques |11[ I12j (see Fig. |^(b)) rather made the assumption that modifying Maxwell’s equa¬ 
tions at one i = iq node would lead to the reflection of the driving wave at this particular node with a 
reflection coefficient r and transmission coefficient 1 — r. This would be equivalent to adding only one pseudo 
source at position with a re-emission coefficient that radiates in both directions i < io and i > Iq. 

While this is adequate for the analysis of second order FDTD schemes for which the former model was 
developed initially HH, the assumption of a single reemission pseudo-source leads to discrepancy between 
modeling and simulations for the coefficient r at orders p > 2 |12) . 



Modified node i 


n 


Figure 2: Principle of the model, (a) The new model takes into account the fact that at order p, several nodes I around 
the modified node act as monochromatic "secondary pseudo-sources" of error with amplitudes ^i, frequency oj (equals to the 
incident wave frequency) and wave vector k given by the dispersion relation in vacuum, (b) Former model considered that only 
the modified node is acting as a monochromatic source of error with amplitude r. 


This inaccurate estimate becomes highly detrimental when looking for instance at stencil truncations 
with the domain decomposition technique at very high orders p [7]. As explained in the above section, the 
discrepancy is due to the fact that for Maxwell solvers of orders p > 2, the 2{p — 2)-adjacent nodes near io 
will also be affected by modification of the fields at ig. As a consequence, they will also act as secondary 
pseudo-sources that radiate monochromatic waves. 

In the following section a method is introduced for computing the reemission coefficients (^ro-i-/)-(p- 2 )sgz^(p- 2 ) 
of secondary pseudo-sources and deduce the total amplitude C of the error generated by their interference. We 
will call this new method the "p-sources" or "multi-sources" model as opposed to the "one-source" models 
developed earlier. 
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2.2. ID Model: Detailed analytical calculation of the coefficients fi 
2.2.1. Initial assumptions 
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Figure 3: Configuration of the problem in a ID geometry. Maxwell’s equations are modified at rim nodes of the grid. As 
shown in the previous subsection, the 2(p — 2) adjacent nodes at left and right of the modified nodes also act as "secondary 
pseudo-sources". 


The problem layout is detailed in Fig. For the sake of simplicity, the method is first detailed in a 
ID geometry for a staggered field solver of order p, and is generalized for a 3D geometry in the following 
subsection. 

We assume a linearly p-polarized {Ed, Bd) driving wave of frequency oj propagating along the x-axis 
(Electric field in the plane (x,y) and Magnetic field in the plane (x,z)). 

We consider an infinite domain along x and assume that Maxwell’s equations are modified at Um + 1 
nodes and n^, internodes of the grid, from i = 0 to z = Um- As a consequence, nodes i and internodes i + \ 
ranging from imin = ~{p ~ to imax = Tim + (p ~ 2) will act as secondary pseudo-sources of amplitudes 
{f,i+l)-{p- 2 )^l^p -2 Each one of this individual source will emit electromagnetic 

radiations in both {x < 0) and {x > 0) directions. 

2.2.2. Expression for magnetic and electric Ey fields 

From the analysis in the preceding subsections, we can infer a form for the total electric and magnetic 
fields {Ef.,B^:) at position i and time n which are the combination of three components for the electric field 
F;, (see Fig. a (a)): 
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Bz,i 
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bbyinternode,i,l-^^ 
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( 21 ) 

( 22 ) 


where the individual components are defined below: 

• the electric field E^^ at position i corresponding to the electric field of the driving plane wave: 

Tpn _ —jckiSx+ujn6t 

■^d,i ~ ^ 


( 23 ) 
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Figure 4: Contribution of adjacent sources and driving wave to the total electric and magnetic fields {Ey, Bz). (a) The expression 
of the electric field Ey^i at node i is a combination of electric field Ed,i at position i of the driving plane wave, the electric field 
CebynodB,i,i ^t position i emitted by the secondary source located at position I {imin < I < imax, the electric field CebyinteTnode,i,i+^ 
at position i emitted by the secondary source located at internode / +1 {imin < Z +1 < imax)- (b) The expression of the magnetic 
field Bz,i at node / is a combination of magnetic field Bd,i at position i of the driving plane wave, the magnetic field ({bbynode,i,i 
at position i emitted by the secondary source located at position I {imin < I < imax, the magnetic field Cbbyintemode i l+^ 
position i emitted by the secondary source located at internode I + | {imin < / + | < imax)- 

• the electric field Qbynode i i position i emitted by the secondary source at position I {imin < I < imax)'- 

n f i ^ I 

^ebynode,i,l | ^^^—jckiSx+ujnSt ^ / (^‘^) 
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(25) 


and similarly three components for the magnetic field (see Fig. 0(b)): 


the magnetic field at position i corresponding to the electric field of the driving plane wave (see 
Fig. Kb)): 

j^n _ ^-jckiSx+umSt ^ 20 ^ 

the magnetic field at position i emitted by the secondary source at internode I + ^ 

{imin <1 + h < 

imax)' 




(27) 
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2.2.3. Calculation of the total error in the general case 

For a staggered field solver of order p, Maxwell’s discrete equations for the electric and magnetic fields 
can be written as: 
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(29) 


and: 
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( 30 ) 
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z,i +5 


“ f^i+h 






9=i '^i+k,q y:*-(<?-i) 


with ai, f3i, and coefficients that may vary with the position on the grid i. In vacuum, we have 
Oi = 1, I3i = —c6t/6x and = V’f; = Cf for all positions i on the grid. 

For the following, let us assume that the coefficients a*, j3i, and are modified between nodes 
i = 0 and i = n^- Under these assumptions, the unknowns of the system of equations with varying 
coefficients are the N = 2{imax — 'imin) + 1 re-emission coefficients coefficients nodes and 

i^l+^)im,iri^l+-^imax mtemodes. Finding them requires the Ne = {imax — imin) + 1 equations (29) for 

min ^i ^ imax and the Nb = {imax - imin) equations (|M1) the 


the electric field Ey written at nodes i 


magnetic field written at internodes imin ^i + 2 ^ imax- Note that the closure of the system is achieved 
thanks to the finite number of reemission sources, as it is assumed that there is no additional sources beyond 
i < imin and i > imax (vacuum). 

In vacuum, it can be shown that the system to be solved can be expressed in the following matrix form: 


with: 


A„X = 0 


(31) 
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(32) 


a band Matrix of size N x N and bandwidth p -|- 1 (total width of the stencil). X is a vector of size x 1 
containing the unknowns of the problem: 


X = 




^i-n 


(33) 


0 and ( 6 ),' . </<,• 
takes the form: 


= 0 . 


Note that equation (31) trivially yields X = 0 (providing that det X 7 ^ 0) in vacuum i.e (^j_|_i )j 


1 ^l+h^i n 


f Maxwell’s equations are modified at Um nodes and Um internodes, the system 
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(35) 
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The linear system of equation (34) can be easily solved analytically using Cramer’s rule. The total 
re-emission coefficient for the electric field in the x > 0 direction is then simply given by: 


^max 

^—'^min 



^—'^min 


(37) 


2.2.4- Analytical resolution of the system 

Analytically, Cholesky LU decomposition m that runs in 0{N^) may be used. However, as the system 
to be solved is sparse, iterative methods for sparse matrixes that scale in 0{NLogN) are used instead. For 
low order systems, the resolution can be done "by-hand" but requires advanced symbolic calculations at 
large orders. Results discussed in the remainder of this paper were obtained using Mathematica [T^ and 
the Linear Solve function that has accelerating techniques for sparse matrixes, allowing fast calculation of 
numerical solutions for the coefficients 
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2.3. Extension of the model to 3D 
2.3.1. Initial assumptions 

We consider a 3D geometry and the case of a driving p-polarized {Ex, Ey, E^, B^) plane wave of 

frequency to propagating with an angle 6 in the (x, y) plane and an angle cj) in the (x, z) plane. In vacuum 
the wave vector k = {kx, ky, kz) = [k cos 9 cos (f, k sin 9, k cos 9 sin (f) of such a wave is simply given by the 3D 
dispersion relation of the staggered scheme. 

In the case of p-polarized driving wave the E field is in the plane of incidence (x,y) and the B field 
is orthogonal to this plane. It can be easily shown that equations are of the same form with a s-polarized 
driving wave where the B field is in the plane of incidence and the E field orthogonal to the plane of 
incidence. Results from arbitrary polarization can then be easily deduced by noticing that any polarization 
is the superposition of these two orthogonal polarization states, justifying the restriction of the analysis to 
the p-polarized case only without loss of generality. 



Figure 5: Problem configuration in the 3D case (i^ = 0°, kz =0). On this illustration, the driving wave is obliquely incident 
with an angle 6 and 0 = 0° on a set of Um modified nodes. Modified nodes act as secondary source planes that radiate in a; < 0 
and a; > 0 directions. 


The configuration of the 3D problem is sketched on Fig. [^in the case (f) = 0° and 9^0. We consider 
an infinite domain along x, y and ^ axes and assume that Maxwell’s equations are modified along direction 
X only (translational invariance along y and z), at Um nodes and Um internodes of the grid, from i = 0 to 
i = nm- In this case, nodes and internodes from imin = ~ 2) to imax = nm + (p ~ 2) will act as a 

secondary pseudo-source planes of amplitudes (?i+/)-(p- 2 )^Zs£(p- 2 ) and (Cj+i+;)_(p_ 2 )^Zsg(p- 2 )- 

As there is no stencil variations along y and z, the spatial phase shift induced by stencil modifications 
will not depend on j and k. This implies that there is no spread in propagation direction of the secondary 
sources. Consequently, each one of this individual source plane will emit plane waves in x > 0 and x < 0 
directions with angles 9 and tt — 9 with regard to the x axis (see Fig. [^. 
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2.3.2. General form of modified Maxwell’s equations in 3D 

In the 3D case, Maxwell’s equations for the electromagnetic fields {E^, Ey, Ez, Bx, Bz) have the following 
general form: 




^x,i+\^x,i+\,j,k 


^^n+ 1/2 




-cB 


n+| 




(38) 


E 


71 +1 


= a. 


rpn 

y,i3+h,k 


+ f^y,i Z)f=i 


cB 


71+1/2 

x,i,j+\,k+\ + {l-l) 


- cB'"^^ 


x,i,j+kk+k-i 


1yd 


T.U ^hcB^ 


p 


«+l/2 _V2 rR 

2 :,i+^+((—l),j+|,fc 


z,i+k-lj+kk 


(39) 


1 

zx,J,k+2 


— I3z,i Z)f=l C'f 


cS 


71+1/2 


x,i,j+h + {l-R,k+h 


1 - 


^dj+h—l-,k+d 


(40) 


cB 


^dj+h,k+k 


* 1 * X,lJ+i fc+i 




+ 


7l,i Sill cf 


- E' 


2:,2,7 +/,/;;+1 2:,2j —(/ —l),fc+|J 

-pn _ pn 

y 4 J+|,fc+^ yd 3 +\,k-{i-i)\ 


(41) 


cB 


z,i+\,j+\,k 


= a 


z,i+ 


1 C-B 


^-2 


24+2>i+2’^ 


- /3 


2 , 1+1 


E 2 T^p pri _ 2 „/,P fpn 

^=1 i+\,l y,i+l,j+\,k ^l=^^i+ki j+l fc 


( 42 ) 


+ ti+iSt.cf 


a;,i+|,j+Z,fc x,i+i,j —(i—l),fc 


o^x,y,z, Px,y,z, Oi^^zi f^x,zGy,zi lx,z coefficients that may vary with index i (not j or k in this configura¬ 
tion). In the above equations the discrete operator is modified by taking different stencil coefficients 
and that may also vary with the position on the grid i. Vy and are not modified along j and k. In 
vacuum, the coefficients are given by: 


(^x,i — *^07,1 — — (^z,i — ^z,i — ^ 

Pz,i = I3*z^i = Jy,i = 7^,1 = Vx 
I3x,i = Px^i = lz,i — 3y 

lx,i = Py,i = 3z 


(43) 

(44) 

(45) 

(46) 

(47) 


and 


7 , = 7 , = Cf ( 48 ) 

In our configuration. Maxwell’s equations between nodes i = 0 and i = Um are modified so that coefficients 
a, (3, 7 , a*, /3*, 7 *, F and if may change between these positions. 
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2.3.3. System of equations to solve in the 3D case 

Due to translational invariance along j and /c, equations (38|40|4f ) can be written as: 


C? _ pjckylSy\ 

pn _ o Z^;=1 '-'I ^ ^ ^r»i+V2 

ik h'x,i+l/2 a ^ 1 h 

2 :,«+ 2 J+ 2 ,/c 

Vf C? _ pjckylSy\ 

TT^n o ^ _ „Dn+l/2 

^ ■ • ' '4 “ “^^4+1/2 j^ojSt - - -- ^ ^ ^ i 1 1 


z,i,j,k+^ 

cB^^.^. 1 ^ 1 


^ 2/,2 + 1/2 




x,i,j+^,k-\--r 


xE^. 1 

2,2 J,/C +2 


C, (7P^p-Jcfc2^52_pJcfc2(^-l)'52'\ 

, * jcUj5t U=i^i{^ _f_ )_ 

^ ^ ^y,i,j+k^k 


(49) 

(50) 

(51) 


By expressing Ex and E^ as a function of Ey and using the above equations (49), (50) and (51), the 


system of equations (38 42) can be re-written into a system of two equations only on Ey and B^'. 


+7' = 


^y,i 


y^2 pP „n"^+l/2 _V2 ih^ rB^^^ 


(52) 


cB'^^K . 1 ^ = a* icS”. / . 1 

£ 


- '9;.,. 


2 pP pn 
^9=1 i+|,g p.*+gj+5 


1 - El., V’L, 


*+5.9 P.*-(9-l).i+5 


with: 


Ui — ^y,i 4” 


+.*7^,4"^"^"^* fell Cf 1 ^ 


gjci^St _ Q, 4 / gici^5t _ Q,* . _ «* gjc{u}St+kySy) y^2 1 Cf (e~Eky(l-^)&y — ejckyldy\ 
yj \ Xf\j ‘ X^ij —1 t \ ) 


1 2 


and: 


oic*^5t 


*^*+5 7^^j_|_i+,j+i _ 




-1 2 


l=l 


(yP ( g-lcfea(^-l)'5p _ gjckylSy 


(53) 


(54) 


(55) 


Equations 52 and 53 having the same form as equations (29) and (30), the method developed with the 


ID model is readily applicable for deriving the re-emission coefficients fi. When 0 = 0 and </> = 0 we have 
ky = 0, kz = 0 yielding exactly the same equations as in the ID case previously described. When 0 7 ^ 0 and 
(j) ^ 0, the form of the total electric field Ey and magnetic field Bz is slightly different in the 3D case and is 
given in the following subsection. 


2.3.4. Form of electric and magnetic fields in the 3D case 

It can be shown using the system of equations ( 38||42 ) in vacuum that the relative amplitudes Xy = Ey/E 
and Xz = BzjB where EjB are the total electric field/magnetic amplitudes are given by: 


2^2 


Xp = 


Xz = 




^(1 - fp)^ + + (1 - e2)2e2e2 


(1 - 


^(l-e 2)2 + e 2 g 2 


(56) 

(57) 

(58) 
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with ex,€y, Cz satisfying the dispersion relation in vacuum: 

and given by: 


_|_ g2 _|_ 2 _ 2 


'^x^y^z 


(21 - 1 ) 


(-X.V.Z — . ^ ^ sin ^ kx.v.z^: 


x,y,z'Jx,y,z 


( 59 ) 


(60) 


sinw5t/2 ^ 

In the ID case {ky = = 0), Cjy = Cz = 0 and Xy = Xz = 1- In the 2D case {kz = 0), = 0, Xy = ^x 

and Xz = 1. 

Retaining the same general form for the total electric and magnetic fields (-E'^j j fc,-B”i j fc) hi the ID 
case: 


i -i 

imax 2 


mn _ mn , 

^y,i,j,k ~ ^d,i,j,k ' 


1 

k,l / ^ebyinternode,iJ,k^l-\-k 

(61) 

1 

^ —'^min 

^ 2 



'^max 

?■ 

•'max 2 


TDU _ on _j_ 

^zd,j,k -^d,i,j,k ' 

^ ^ Cbbynode,iJ,, 

k-,^ / ^bbyinternode^iJ,k,l^^ 

(62) 


^—'^min 

^—'^min 2 

(63) 


the combination of three components for the electric field Ey (see Fig. (a)) is now given by: 

• the electric field E^^j ^ at position {i,j, k) corresponding to the electric field of the driving plane wave: 




_ —jckxiSx—jckyj6y—jckzk5z-\-jci^n5t 


(64) 


the electric field Qbynodeijkl position {i,j,k) emitted by the secondary source plane at position I 

{‘^min < i < imax)' 


Qbynode,i,j,k,l 


-^^^^g^jckxiSx-jckyjSy-jckzkSz+jcionSt 

^^^^^-jckxiSx-jckyidy-jckzkSz+jcLonSt 


i ^ I 
i > I 


(65) 


• the electric field ^ ^ i nt position (i,j,k) emitted by the secondary source plane at in¬ 
ter node I 2 iXmin ^ 2 ^ ^max) • 


c 

^ebyinternode,iJ^k,l-\- 


1 

2 


~Xy^l -\-1 

Xy^l 1 Q-jckxiSx-jckyi5y-jck^k6z+jcLonSt 


i ^ I 
i ^ I 


while the three components for the magnetic field Bz read: 


( 66 ) 


• the magnetic field at position {i,j,k) corresponding to the electric field of the driving plane 

wave: 


on _ -jckxiSx-jckyjSy-jckzkSz+jcOJnSt 
^d,i,j,k - XzC 


(67) 


the magnetic field C 


1 + ^ {i min <i + l < 

imax)‘ 


bby internode,i,j,k,l+k 


at position (i, j, k) emitted by the secondary source at internode 


^bhy internode,i,j,k,l+k 


Xz^l 1 gickxiSx-jckyjSy-jckzkSz+jcixinSt 


Tz .f 1 p-jckxiSx-jckyjSy-jck^kSz+jcUmSt 


i ^ I 
i ^ I 


( 68 ) 


15 




the magnetic field Cbbynode,i,j,k,i Position 


{^min < i < imax)' 


{i,j,k) emitted by the secondary source at internode I 


Cbbynode,i,l 


^^^j^^-jckxiSx-jckyj5y-jckzk5z+jcU)n5t 


i ^ I 
i ^ I 


(69) 


3. Numerical validation 1: application of the model to the domain decomposition technique 

We now investigate the effect of domain decomposition on solutions calculated with high-order/pseudo- 
spectral solvers. 

3.1. General principle of the domain decomposition technique 

The principle of the domain decomposition technique is illustrated on Fig. For simplicity, we consider 
first a ID geometry along x. The main domain is split into two semi-infinite subdomains (Di) for x < 0 
and {D 2 ) for X > 0 with Ng^ards guard cells (Gi) = [0, Ng^ards^x] for (Dl) and Ng^ards guard cells (G 2 ) = 
[—guards^x^D] for {D 2 ). At each time step n, the procedure for exchanging data between the two domains 
is as follow: 



Figure 6: Domain decomposition technique. The simulation domain is divided in two semi-infinite subdomains (Di) and (D 2 ). 
The envelope of the stencil of the field solver (order p) is represented by a black curve and grey filling between the axis and the 
curve. At each time step n, fields are first updated using Maxwell’s equations on each subdomain. Fields in guard regions (Gi) 
and (G 2 ) are then replaced by fields at the same positions in (D 2 ) and (Di) respectively (red arrows). 


(1) the electromagnetic field component is updated using Maxwell’s equations independently on each sub- 
domain (Dl) and (D 2 ) (Fields in (Gi) and (G 2 ) are not updated), 

(ii) nodes from (D 2 ) at same positions as nodes in (Gi) are copied to (Gi), 

(hi) nodes from (Di) at same positions as nodes in (G 2 ) are copied to (G 2 ). 

This procedure is done for every field component. Using finite-difference leapfrog time integration, 
is first advanced to and copied to adjacent guard regions. is then advanced to using 

and copied to adjacent guard regions. If Nguards ^ p/2, domain decomposition will not affect the precision 
of the scheme and the result will be equal to the one on a single grid without domain decomposition. 

Because exchanging large volumes of data can be prohibitively expensive in terms of computational 
ressources on massively parallel supercomputers, it has been proposed to limit the volume of data exchanges 
between subdomains, and thus the number of guard cells for high-order stencil, or its infinite order pseudo- 
spectral limit |7]- However, errors are introduced when Nguards < p/2. In this case, pj^l — Nguards nodes 
are affected in (Di) and (D 2 ) by truncation of the stencil. This eventually leads to spurious signals in each 
subdomain and potential numerical errors. In the following we use our model to compute the total error 
generated by field truncations and exchange at the boundaries and compare theoretical results to simulations 
for plane waves at normal incidence or at an angle from the domains’ interface. 
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3.2. Domain decomposition strategies for staggered field solvers 

In the numerical validations of our model, we will compare two domain decomposition techniques for 
staggered field solvers. The first one (staggered method) is sketched on Fig. [^(a) for a ID geometry. E and 
B fields are exchanged in a staggered way. 
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(D) 

1 



(a) 

(D ) 
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i = 

=0 

(b) 

(D ) 

2 

E 

1 1 

: i 

(D ) 

2 


B 



B _i _ ; 

I t 


i = l/2 


i = l/2 


Figure 7: Domain decomposition strategies, (a) First domain decomposition strategy (Strategy 1). E and B fields are exchanged 
in a staggered way. (b) Second domain decomposition strategy (Strategy 2) where E and B are exchanged in a centered way. 
The value of the electric field at the central node i = 0 is averaged by taking Efff = + Efff jjfi/2, where 

(resp. is the E-field value calculated on Di (resp. D 2 ) at i = 0. 


The second technique (centered method) consists of centering the exchange of E and B fields by overlap¬ 
ping nodes from {Di) and (D 2 ) (see Fig. Id (b)). In that case, the value of the electric field at the central node 
i = 0 is averaged by taking where (resp. is the E-field value 

calculated on Di (resp. D 2 ) at z = 0. In the following we will refer to the staggered method as "Strategy 
-1" and the centered method as "Strategy-2". 


3.3. Practical implementation of the domain decomposition technique in our model 

For "Strategy-1", the domain decomposition technique is accounted for using the following coefficients 
Maxwell’s equations (29 30) for the ID model and (52 |^) for the 3D model: 


m 


pP - 




rf, = 


Ci,p 

i > 0 .or. p/2 < Nguards 


Ci,p 

i <C 0 .and. i 1 ^ ^guards 


0 . 

i <C 0 .and. i 1 ^ ^guards 

(70) 

Ci,p 

i — 0 .and. -t- / ^ ^guards 


0 . 

i — 0 .and. i 1 ^ ^guards 


Ci,p 

i ^ 0 .or. p/2 <C 



i ^ 0 .and. % 1 ^ ^guards 

(71) 

0 . 

i 0 .and. i 1 ^guards 


rosition technique is accounted for using the following coefficients: 

Ci,p 

Z ]> 0 .or. p/2 <C ^guards 


Ci,p 

i <C 0 .and. i -1- / ^ ^guards 


0 . 

i <C 0 .and. i 1 ^ ^guards 

(72) 

Ci,pl2 

i 0 .and. Z -1- / ^ ^guards 


0 . 

% — 0 .and. % 1 ^guards 
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Besides, in both cases we have: 


Cl^p i < 0 .or. 'pI‘1 < ^guards 
^l,p i ^ 0 .and. i I ^guards 

0 . i>0 .and. i-l < Ng^ards 

Cl^pl‘^ i = 0 .and. i / ^ ^guards 

0 . i 0 .and. i I ^guards 


and 


a 


xg, 


= a 


X,l 


— — ^z.i — 


z^i 


* 2,2 


= 1 


( 73 ) 


(74) 


(^z,i ~ l^z,i — ly^i — ly,i — Vx 

Px,i = = ll,i = Vy 

7x,i = f3y,i = Vz 


(75) 

(76) 

(77) 

(78) 


3.4- Theory-simulation comparisons 

In this section, we validate our model against simulations in the case of a plane monochromatic wave of 
frequency oj impinging at oblique incidence {9 ^ 0°,4> = 0°) on the subdomain boundary. 

3.4- 1- Test procedure and numerical error measurements 

The simulation domain is divided into two subdomains of equal sizes with Ng^ards cells at their boundaries. 
At each time step, fields are computed independently on each subdomain and guard cells are exchanged 
between subdomains. At t = 0, a Harris-like waveform H{t,y) is launched at the left boundary of the 
simulation domain: 


H(t,y,X) = h{t) sin {cot — ky sin 6) (79) 

where 6 is the angle of incidence of the waveform, k =‘lie jX the wavenumber and uo the angular frequency 
obtained by the numerical dispersion of the staggered scheme. Amplitude /i(t) is the Harris function given 
by: 


h{t) 


^ [lO — 15 cos ^ 6 cos ^ — cos 

0 


0 < t ^ T 
t > T 


(80) 


This wave-form has a quasi-monochromatic spectrum, enabling validation of the model at specific wave¬ 
length A. The modulus of the numerical re-emission coefficient |Csim| is obtained by taking the ratio of the 
reflected field energy R over the incident wave energy I: 


,, , ,[R jELeftiE^+C^B^) 

~ y j ~ y j 


(81) 


where the sum ^ in the above equation is taken on the left subdomain and at the end of the simulation. 
In order to obtain also the phase properties of the secondary sources in the simulations, we compute |1 —CsimI) 
which equates the ratio of transmitted energy T over the incident wave energy T. 


|1 Csiml — \ T — 



EmgHtiE^ + c^B^) 


(82) 


where the sum ^ in the above equation is taken on the right subdomain and at the end of the simulation. 
The phases Arg{^sim) and Arg{l — (sim) of the reflected and transmitted waves in these simulation are Anally 
given by: 
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Arg{Csim) 


Argi^l Csim') 


arccos 


arccos 


1 + ICsimP ~ |1 ~ G 
2|Gim| 

1 + 11 ~ Gimp ~ IG 
2|1 - GimI 


(83) 

(84) 

(85) 


In the following paragraphs, we compare, for different numerical parameters, the coefficients |Gim|) 
|1 - GimI, Arg{Csim) and Arg{l - C,sim) to our theoretical estimates \Cth\, |1 - Cth\, Arg{Cth) and Arg{l - C,th) 
provided by the model described in section 

3.4-^- Influence of order p and number of guard cells Ng^ards 



Figure 8: Variation of the truncation error C, (total re-emission coefficient) with order p for a different number of guard cells 
Nguards = 3,5,10,20,30 corresponding to red, blue, green, magenta and black curves. X/5x = 10, and {9 = 0°,(j> = 0°) are 
fixed. "Strategy-2" was used as domain decomposition technique. For each curve, circles represent the re-emission coefficient 
calculated in the simulations, solid lines represent ^ computed with our model and dotted lines represent |(C| deduced from the 
analytical approximate formula derived in [Appendix B.lj 


Figure]^ (a) represents the variation of the modulus of the total re-emission coefficient |G with order 
p of the Maxwell solver for different number of guard cells Nguards (see colored curved). The agreement 
between the model (solid lines) and simulations (circles) is excellent. When Nguards ^ there is no 
stencil truncation and |^| is zero at double machine precision Xdouble ~ 10“^^. However, when Nguards < p/2 
the stencil is truncated at subdomain boundaries and a spurious signal of amplitude |((| is created. 

For a given number of guard cells Nguards, the stencil truncations increase when order p increases and the 
amplitude f also becomes larger. Figure]^ shows that seems to converge to a a constant value when p S> 1 
allowing the estimation of Coo when p —oo i.e when the FDTD Maxwell solver turns into a pseudo-spectral 
solver. 

On the contrary, when Nguards increases for a given oder p, stencil truncations are reduced and |C| 
decreases, as expected. 

Note that it is possible to get a practical analytical approximate formula of the error amplitude |C| (cf. 
Appendix B.l). This estimate is represented by dotted lines on Fig. [^and dashed lines on Fig. It very 


accurately reproduces the evolution of |C| with p, Nguards 5x/X. 

The limit |C|oo of this analytical approximation when p —>■ oo is derived in Appendix B.2 The variation 
of the error |C|oo with the number of guard cells Nguards is represented in Fig. (a) and (b) (green dashed 
line). The maximum error (green line) varies as when Nguards 1 which could have been guessed 

by simply noticing that the amplitude of the stencil coefficients IC^j vary as 1/(2/ — 1)^. This means that 
taking 10 times more guard cells yields a maximum error divided by 100. Even if the absolute error is 
relatively low, having zero error at machine precision (~ 10“^^ for double precision) would mandate too 
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Figure 9: Variation of the truncation error (total re-emission coefficient) with the number of guard cells for different orders 
p and wavelengths X/Sx. "Strategy-2" was used as domain decomposition technique, rjx = 0.4 and {6 = 0°,(p = 0°) are 
fixed. For each curve, dashed lines represent the re-emission coefficient |((| deduced from the analytical approximate formula 
(cf. [Appendix B.l l. Red, Blue and Magenta circles represent |^| computed with the full analytical model. The green solid 
line represents the function (a) Variation of |C| with Nguards for \/5x = 5 (b) Variation of |^| with Nguards for 

X/5x = 100. 


many guard cells. Fig. shows however that for orders as high as p = 1000, the number of guard cells 
needed to have zero error at machine precision remains very low {Nguards ~ 100 for p = 1000). As order 
p = 1000 already yields spectral resolution at machine precision on a large band of frequencies, it would 
thus be more interesting to use very high finite order solvers instead of infinite order solvers with domain 
decomposition. These results will be presented in greater details in further work. 


3.4-3. Influence of wavelength A of the driving wave 

The model is now compared to simulations when X/5x is varied for fixed parameters Nguards, Vx = 0.4 
and 9 = (j) = 0°. Results of the model and simulations are represented in Fig. 10 Our analytical model 


again perfectly reproduces simulation results. The analytical approximate {x markers in Fig. 10 (a)) derived 


m 


Appendix B.l again yields very accurate estimates of the error magnitude. 


The variation of the modulus |((| and phase Arg{f) of the total re-emission coefficient are represented on 
panels (a) and (b). As expected decreases for longer wavelength. This can be qualitatively understood by 
considering the fact that for long wavelengths A, the stencil is truncated on a region that becomes spatially 
very small compared to A. For any wavelength, the emitted field in the x < 0 direction is dephased by a 
constant tt/2 from the driving field, which is actually equivalent to a reflection of the incident field on the 
subdomain boundary with a reflection coefficient (4- 

The modulus |1 — Cl and phase Arg{l — C) of the transmitted wave in the x > 0 direction are represented 
on panels (c) and (d). Panel (c) shows that the amplitude |1 — Cl of this wave is approximately equals to the 
amplitude of the driving wave on the whole frequency domain and with a very low dephasing. Consequently, 
there is low effect of stencil truncations on the wave passing through subdomain boundaries. 

However, our model shows that energy created in the simulation box AE after boundary crossing is given 
by: 


AEfEinc 


(|1-CI + ICI)-1 

(86) 

|l_|C|ef-^"9(0| + |c| -1 

(87) 

|i-jc|CII + ICI -1 

(88) 

vFW+ICI-i 

(89) 

i + 1T + |c|-i + „(|c|2) 

(90) 

ICI + o(ICI) 

(91) 
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Figure 10: Variation of the truncation signal amplitudes and phases with Sx/X. Parameters Ng = 5, = 0.4 and 9 = (j> = 0° 

are fixed. "Strategy-2" was used as domain decomposition technique. The blue lines correspond to the case p = 20 and the 
red lines to p = 12. In each plot, solid lines correspond to results of the model. " + " and "o" markers correspond to simulation 
results, "x" markers correspond to deduced from the analytical approximate formula (cf. Appendix B.ll (a) Variation of 
the amplitude of the reflected signal as a function of the ratio X/Sx. (b) Variation of the phase Arg{^) of the reflected signal 
as a function of the ratio X/Sx. The dashed lines correspond to relative errors in % between the model and simulations for 
p = 20 (blue lines) and p = 8 (red line), (c) Variation of the amplitude |1 — Cl of fh® transmitted signal as a function of X/5x. 
The dashed lines correspond to relative errors in % between the model and simulations for p = 20 (blue lines) and p = 8 (red 
line), (d) Variation of the phase Arg{l — C) of the transmitted signal as a function of X/Sx. 


10 (a)), Arg{C/) k, 7r/2 (cf. Fig. ^(b)) and Einc the initial energy of 


where we used <C 1 (cf. Fig. 
the driving wave. 

The total energy after boundary crossing is thus now greater than the initial energy of the driving wave 
and equals to |C|- This means that stencil truncations created non-physical energy in the simulation box of 


magnitude |C| (cf. Fig. 10 (a)). If the driving pulse crosses several boundaries during the simulation, the 


total energy will thus increase through time with a growth rate that is a function of the number of crossed 
subdomains Ngub |(C|. Moreover, if the number of crossed subdomains is high so that Ngui,Arg{\ — ~ 1 

dephasing effects could also start to alter the phase of the transmitted wave. These effects are under study 
and further analyses will be presented in greater details elsewhere. 

3.4-4- Influence of the angle of incidence 9 

Fig. 11 shows the variation of the modulus of the spurious reflected signal (relative to the incident 


wave modulus) |C| with the angle of incidence 6 for two different wavelengths \/6x = 3.4 (red curves) and 
A/(5x = 10 (blue curves). In all cases, it appears that variations of |C| with 9 are small for low angles. 

Notice that in our case, there is no stencil variations along directions y (index j) and z (index k). This 
implies that the spatial phase shift induced by the domain decomposition is independent of j and k. As 
a consequence, there is no spread in propagation direction of the reflected/transmitted waves, which is in 
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Figure 11: Variation of the truncation signal amplitude (re-emission coefficient) |i(| with the angle of incidence 6 for the two 
different X/5x = 3.4 (red curves) and \/5x = 10 (blue curves). Parameters Ng = 5, rjx = 0.4 and order p = 20. "Strategy-1" 
was used as domain decomposition technique. Solid lines represent results of the model and points results from simulations. 


accordance with the initial assumptions of the 3D model. 
3.4-5. Influence of the domain decomposition strategy 
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A/6x 


Figure 12: Variation of the truncation signal amplitude Id with 5x/X for the two different domain decomposition strategies. 
Parameters Ng = 5, px = 0.4, order p = 20 and 0 = (j> = 0° are fixed. Red curves correspond to "Strategy-1" and blue curves 
to "Strategy-2". Solid lines represent results of the model and points results from simulations. 


Fig. 


12 illustrates the effect of the domain decomposition method on the re-emission coefficient |(^|. 

1^1 remains fairly uniform on the whole spectral range. On the 
I Cl decreases with X/6x and is signihcantly lower than 
oo. This indicates that Strategy-1 method 


In the case of "Strategy-1" (red curves) 
contrary, in the case of "Strategy-2" (blue curves 
with " Strategy-1" at large X/6x and seems to vanish as A 
induces a systematic error that is not present with Strategy-2 and that the latter is thus preferable. 


4. Numerical validation 2: application of the model to Perfectly Matched Layers 

In this section we illustrate the versatility of the new analytical model through its application to the 
evaluation of the reemission coefficient |C| of a PML. Former studies noticed a discrepancy between theory 
and simulations using a ’one-source’ model at high orders p [12]. In the following, simulations and results of 
this ’one-source’ model are compared to om "multi-sources" model. 

4-1. Practical implementation of a PML in our model 

Numerical parameters for the PML are the same as in |12j : 
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• electric/magnetic conductivities in the PML of a(i} = amaxii^x/A)'^ with amax = 4/(5x and A = 55x, 

• the PML is made of Np^i = 20 cells. 

As shown in Appendix B, a PML (using Berenger’s original split formulation in this example) is described 
in our model by using the following coefficients (see Appendix B for detailed calculations): 


with: 


as well as: 


€i = rf, = Q,pyi 


0'y,i = 


Ct . , 1 


0^z,i — (^x,i ~ 1 

a{i) 

1 


l^x,i+^ 

* 

To;,2 

I3y,i 




(92) 


(93) 

(94) 

(95) 

(96) 


= 

= ^*x,i = = Vy 

= Vz 

= m 


(97) 

(98) 

(99) 

( 100 ) 

( 101 ) 

( 102 ) 


and 


a{i) 

P{i) 

7(i) 


2 — a{i)5t 
2 + a{i)6t 
2c2 


iVx 


2 + a{i)6t 

^jcuiStl2 _ ^-jciLiStf2 

— Q:(z)e 


(103) 

(104) 

(105) 


In the following, our "multi-sources" model is used to compute the total re-emission coefficient from a 
PML and is compared to theoretical and simulation results of [12] that uses the "single-source" model. 

4-2. Influence of wavelength A on 


Fig. 13, represents the variation of |C| with the wavelength X/6x. Panels (a) and (b) of Fig. 13 show 
that the new "multi-sources" model agrees perfectly with simulation results, validating its predictive value. 
Furthermore, it explains and solves the origin of the discrepancy observed between the one-source model and 
simulations at high orders. Indeed, Panel (b) shows a disagreemen10 of about 20% between the "1-source" 
model and the simulations results at order p = 8 (blue dashed line). As explained before, this mismatch 
between theory and simulation comes from the fact that previous models only considered one secondary 


^Notice that the discrepancy between the "1-source" model and simulations seem lower than 20% on panel (a) because data 
are plotted in log-scale. 
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source at high orders. Our model however solves this discrepancy (red line on panel (b)) by including all the 
secondary sources that also radiate and couple each other to contribute to the total re-emission coefficient 

IC|. 




Figure 13: Total re-emission coefficient of a PML (a) Comparison of the total re-emission coefficient computed with 
our analytical model (red line) and with ID numerical simulations (red dots). Fixed parameters are the courant parameter 
77 = St/Sx = 0.4, angle of incidence {9 = 0,(j) = 0) and the order of the field solver p = 8. The dashed blue line represents 
the reflection coefficient computed with the "1-source" model developed in |12| (b) The ref line correspond to the relative error 
1C ~ Csimu\/\Csimu\ in % between C calculated by our "p-source" model and (simu calculated by numerical simulations. The red 
line correspond to the relative error |C — Caimul/Csimu in % between C calculated by the " 1 -source" model as developed in m 
and C,simu calculated by numerical simulations. 


4-3. Influence of angle of incidence 9 of the driving wave on the PML 



K-T 



0 (rad) 


Figure 14: Variation of the total re-emission coefficient of a PML |C| with the angle of incidence 9. Fixed parameters are 
g = 6t/5x = 0.4, (j) = 0 and order p = 64 (a) Variation of |C| with 9 (in rad) for \/Sx = 4. The red line corresponds to 
Id computed with our "multi-sources" model and the red points to |Csim| measured in simulations. The blue dashed line 
corresponds to |C| computed with the "1-source" model as developed in | 12 | . (b) Variation of |C| with 9 (in rad) for X/Sx = 8 . 
The red line corresponds to |C| computed with our "multi-sources" model and the red points to |Csim| measured in simulations. 
The blue dashed line corresponds to |d computed with the "1-source" model as developed in |12| . 


Fig. represents the variation of the re-emission coefficient |(^| with the angle 9 for two different 
wavelengths \/6x = 4 (Panel (a)) and X/6x = 8 (Panel (b)). The solid red line corresponds to our "multi¬ 
sources" model and the dashed line to the "1-source" model. In both cases, our model perfectly matches 
simulation results. 
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5. Conclusion and prospects 

This paper presented a novel very general approach allowing the prediction of errors induced by any mod- 
ihcations of stencil in discrete solvers of Maxwell’s equations on cartesian grids. The scope of this method is 
broad as it in principle applies to any linear system of discrete equations. In the case of electromagnetic sim¬ 
ulations, our study demonstrates that this model can be efficiently used to predict the amount of truncation 
errors coming from the use of domain decomposition technique or PML with very high-order/pseudo-spectral 
solvers. 

In particular, our model shows that very high order held solvers can still be used with the cartesian 
domain decomposition technique and a reasonably low number of guard cells Nguards without signihcant loss 
of accuracy. The number of guard cells Nguards so that truncation errors Q do not spoil the required accuracy 
in the simulation can now be predicted accurately with the new model. 

Besides, the general formalism provided by this approach allows the implementation of any conhguration 
and the testing of arbitrary domain decomposition technique at arbitrary order. For instance, this model 
may be used before running a large 3D parallel electromagnetic simulation to predict the total amount 
of truncation errors coming from the use of domain decomposition/PML and choose the optimal set of 
parameters to obtain the hnal solution with a given accuracy. This is of great interest as running our model 
will be vastly faster than running the actual 3D simulation to effectively measure truncation errors and adapt 
the numerical parameters afterwards. 
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Appendix A. Prom arbitrary-order solvers to pseudo-spectral solvers 


We briefly review the arbitrary-order staggered Maxwell’s scheme which is amongst the most commonly 
used solvers in electromagnetic simulations. 

We then verify analytically that when the order p of this scheme tends to infinity, the arbitrary order 
solver converges to a pseudo-spectral solver jlOj . 


Appendix A.l. Arbitrary order seheme 

Discretized Maxwell’s equations in space and time in vacuum can be written in the following general 
form: 

= E” - c<5tv; X (A.l) 

^ cBr^-c<5tV; X E(f (A.2) 

where n is the time step, E the electric field and B the magnetic field, r is the discrete grid on which 
components of the electric field E are defined, s is the discrete grid on which components of the magnetic field 
B are defined and V* is the discrete finite difference operator of order p. For centered schemes components 
of E and B are defined at the same grid points but are shifted spatially for staggered schemes. 

For instance, if we take the simple ID case of a linearly polarized (Ey,Bz) wave propagating along the 
X — axis, we get the following equations for an arbitrary order staggered scheme of order p: 


TTin+l 


2,1+2 


^;‘..-wEcr 


1=1 


^"■+2 

z,i+\ + {l-l) z,i+\-l 


cB^-^ 

2,2 + 


c5t 
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Tpn rpn 


i=l 


(A.3) 

(A.4) 


where i is the discrete position along the x-axis, 5x the mesh size and Cf the coefficient of the finite 
difference operator V* for staggered schemes. These coefficients were first heuristically derived by Fonberg 
|15| . Using the approach developed in usmz], it is possible to derive a closed-form for the Cf coefficients in 
the particular case of the staggered scheme: 


When p —?• oo this yields: 


Cf 


(_l)Uii6i 2(p—l)!^ 

( 2 Z _ 1 ) 2 ( 1 +/_ 1 )!(|_/)!( 2 _ 1),2 
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Appendix A.2. Infinite limit: pseudo-spectral solver 

When the order p —)■ oo, the arbitrary order solvers converges to a pseudo-spectral solver. Fourier 


transforming equations (A.l A.2) with respect to spatial coordinates yields: 



= e'" - c5tV;.cB”^^ 

(A.7) 


= cB"“^- c5tV;.E” 

(A.8) 



(A.9) 
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with E, B and V* the Fourier transform of E, B and V*. Notice that the spatial staggering of elec¬ 
tromagnetic components is implicitly included in the expression of V*. Each component is 

simply given by: 




1 


(1 - ls)26m ^ 






(A.IO) 


where Ig = 1 /2 if the components of the electromagnetic fields in the direction of the finite difference are 
staggered (/^ = 0 otherwise) and = —1. When p —)• oo, we find: 




E 

In the case of the staggered scheme, we get: 




<»(*».) = -;§-E 


(- 1 )^ . km6m{2l - 1 ) 


sm ■ 


7r6m ^ (21 - 1)2 2 

which is the Fourier series development of a triangle wave. For \km^rn\ < this yields: 

V;o(k)=jck 

and Maxwell’s equations write: 

^n+l - _ ^n+1/2 

E = E -jckc^tB ' 

^ n +1/2 1/2 _ _ 

cB = cB — jckc(5t.E 


(A.ll) 


(A.12) 


(A.13) 


(A.14) 

(A.15) 

(A.16) 


which are the equations of the PSTD staggered scheme |10) . Our study at arbitrary order p is thus very 
general as it will also give us information on the behavior of pseudo-spectral solvers when p —)• oo in presence 
of stencil modifications. 
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Appendix B. Analytical approximate of the total re-emission coefficient ^ 

Here we provide an accurate analytical approximate of the total re-emission coefficient ( for 0 = 0 and 
^ = 0 (ID case). The approach can be easily generalized to the 3D case (not detailed here). 


Appendix B.l. Finite order p 

By noticing that stencil coefficients Cf decrease fast with I, we can infer that the pseudo-sources ampli¬ 
tudes due to stencil truncations will also decrease fast from the subdomain boundary where the maximum 
truncation occurs. Considering only a finite number uq < 2{p — 2) of pseudo-sources (instead of the required 
2{p — 2) pseudo-sources at order p to get the exact analytical solution) in the calculation of Q should thus 
yield an accurate analytical approximate of Q. Practically, we checked that using = 3 pseudo-sources .^_i, 


^0 and ^1 in the calculation of C (cf. equation (37)) yields accurate estimate of C- 
2 


C = + ^0 - 
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where ^_i, and are solution of equation (34) with Am = and B = ( 6 i)i<gi<g 3 given 
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and: 
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(B.ll) 
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where "Strategy-2" was used as domain decomposition strategy. 
Solving equation (34) yields: 
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Replacing expressions of C_i) Co and (^1 given by equations (B.19), (B.20) and (B.21) in equation (B.l) 


gives an accurate analytical approximate of C- 


Appendix B.2. Infinite order p —)• 00 
When p —7- 00 , we have: 
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and the elements of matrix Am and vector B can be written as: 


(B.18) 
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TT 


where: 

(i) ^l(^, s, a) is the Lerch transcendent dehned as ^l{z, s, a) = YlV=o + a)“*. This function can be 
evaluated to arbitrary numerical precision in Mathematica, 

(ii) pFq{a,b, z) is the generalized hypergeometric function which can be evaluated to arbitrary numerical 
precision in Mathematica. 


Solving equation (34) at inhnite order then gives: 
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and C,oo can then be calculated using: 
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Appendix C. PML modeling with our multi-sources model. 

The standard formulation of the PML [18] for a p-polarized wave (with 9^0 and (f> = 0) requires the 
splitting of Bz on two components B^x and Bzy that verify: 
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with Bz = Bzx + Bzy and a, (3 coefficients of the PML. In the case of a plane monochromatic driving 
wave of frequency ui, replacing the expression of Bzx and Bzy in Bz we get the following equation for Bz'- 
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which is equivalent to take the following coefficients in the general formulation of equation (42): 
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